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Abstract. We find the numbers of 3 x 3 magic, semimagic, and magilatin squares, as 
functions either of the magic sum or of an upper bound on the entries in the square. Our 
results on magic and semimagic squares differ from previous ones in that we require the 
entries in the square to be distinct from each other and we derive our results not by ad hoc 
reasoning but from the general geometric and algebraic method of our paper "An enumer- 
ative geometry for magic and magilatin labellings" . Here we illustrate that method with a 
detailed analysis of 3 x 3 squares. 
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1. Introduction 

"Today, the study of magic squares is not regarded as a subject of mathematics, but many 
earlier mathematicians in China and Japan studied it." These words from Shigeru's history of 
old Japanese mathematics [TU p. 435] are no longer completely true. While the construction 
of magic squares remains for the most part recreational, their counting has become part of 
the mainstream of enumerative combinatorics, as an example of quasipolynomial counting 
formulas and as an application of Ehrhart's theory of lattice points in polytopes. There are 
several classical [121 El [IS] cind recent [H |2] mathematical works on counting something like 
magic squares, but without the requirement that the entries be distinct, and often omitting 
the diagonals. In previous articles [5l |6] we established the groundwork for an enumerative 
theory of magic squares with distinct entries. Here we apply those geometrical and algebraic 
methods to solve the problem of counting three kinds of magical 3x3 squares. 

Each square x = {xjk)3x3 has positive integral entries that satisfy certain line-sum equa- 
tions and distinctness conditions. In a weakly semimagic square every row and column sum 
is the same (their common value is called the magic sum); in a weakly magic square each of 
the two diagonals also adds up to the magic sum. Such squares have been studied before (see, 
e.g.. Beck et al. [2] and Stanley [IZ]); the difference here is that we count strongly magic or 
semimagic squares, where all entries of the square are distinct. (Since strongly magic squares 
are closest to what are classically known as "magic squares" — see the introduction to our 
general magic article [B] — we call strong squares simply "magic" or "semimagic" without 
qualification.) The third type we count is a magilatin square; this is a weakly semimagic 
square with the restriction that the entries be distinct within a row or column. The numbers 
of standard magic squares (with entries 1,2,..., n"^) and latin squares (in which each row or 
column has entries 1, 2, . . . , n) are special evaluations of our counting functions. 
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We count the squares in two ways: by magic sum (an affine count), and by an upper 
bound on the numbers in the square (a cubic count). Letting N{t) denote the number of 
squares in terms of a parameter t which is either the magic sum or a strict upper bound on 
the entries, we know by our previous work [6j that is a quasipolynomial, that is, there are 
a positive integer p and polynomials Nq, Ni, . . . , Np^i so that 

The minimal such p is the period of N; the polynomials A^^o, Ni, . . . , Np_i are the constituents 
of N, and Nq is the principal constituent. Here we find an explicit list of constituents 
and also the explicit rational generating function N(x) = X]t>o ^{t)^^ (from which the 
quasipolynomial is easily extracted). 

Each magic and semimagic square also has an order type, which is the arrangement of 
the cells in order of increasing value of their entries. The order type is a linear ordering of 
the cells because all entries are distinct in these squares. There are 9! = 362880 possible 
linear orderings but only a handful are order types of squares. Our approach finds the actual 
number of order types for each kind of square; it is the absolute value of the constant term 
of the principal constituent, that is, |A'"o(0)|. (See Theorems 3.4 and 3.14 and Examples 
3.11, 3.12, and 3.21 in our paper on magic labellings [6].) Obviously, this number will be the 
same for cubic and affine counts of the same kind of square. (There is also an order type for 
magilatin squares, which is a linear ordering only within each row and column. As it is not 
a simple permutation of the cells, we shall not discuss it any further.) 

One of our purposes is to illustrate the technique of our general treatment ||6j . Another is 
to provide data for the further study of magic squares and their relatives; to this end we list 
the exact numbers of each type for small values of the parameter and also the numbers of 
symmetry types, reduced squares, and reduced symmetry types of each type (and we refer to 
the On-Line Encyclopedia of Integer Sequences (OEIS) [15j for the first 10,000 values of each 
counting sequence). A square is reduced by subtracting the smallest entry from all entries; 
thus, the smallest entry in a reduced square is 0. A square is normalized by being put into 
a form that is unique in each symmetry class. Clearly, the number of normalized squares, 
i.e., of equivalence classes under symmetry, is fundamental; and the number of reduced, 
normalized squares is more fundamental yet. 

There are other ways to find exact formulas. Xin [IB] tackles 3x3 magic squares, counted 
by by magic sum, using MacMahon's partition calculus. He gets a generating function that 
agrees with ours (thereby confirming both). Stanley's idea of Mobius inversion over the 
partition lattice [ITi Exercise 4.10] is similar to ours in spirit, but it is less fiexible and 
requires more computation. Beck and van Herick [3] have counted 4x4 magic squares using 
the same basic geometrical setup as ours but with a more direct counting method. 

Our paper is organized as follows. Section |5] gives an outline of our theoretical and com- 
putational setup, as well as some comments on checks and feasibility. In Section E] we give 
a detailed analysis of our computations for counting magic 3x3 squares. Sections S] and O 
contain the setup and the results of similar computations for 3x3 semimagic and magilatin 
squares. We conclude in Section [6] with some questions and conjectures. 

We hope that these results, and still more the method, will interest both magic squares 
enthusiasts and mathematicians. 

2. The technique 
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2.1. General method. The means by which we solve the specific examples of 3 x 3 magic, 
semimagic, and magilatin squares is inside-out Ehrhart theory [5J. That means counting the 
number of fractional points in the interior of a convex polytope P that do not lie in any 
of a certain set "K of hyperplanes. The number of such points is a quasipolynomial function 
Epar^it), the open Ehrhart quasipolynomial of the open inside-out polytope {P°,'K). The 
exact polytope and hyperplanes depend on which of the six problems it is, but we can describe 
the general picture. First, there are the equations of magic; they determine a subspace s of 
all 3 X 3 real matrices which we like to call the magic subspace — though mostly we work in 
a smaller overall space M.'^ that results from various reductions. Then there is the polytope 
P of constraints, which is the intersection with s of either a hypercube [0, 1]^ or a standard 
simplex {x G : x > 0, Xl^ifc ~ -'-}• former when we impose an upper bound on the 
magic square entries and the latter when we predetermine the magic sum. The parameter t 
is the strict upper bound in the former case (which we call cubical due to the shape of P), the 
magic sum in the latter (which we call affine as P lies in a proper affine subspace). Finally 
there are the strong magical exclusions, the hyperplanes that must be avoided in order to 
ensure the entries are distinct — or in the magilatin examples, as distinct as they ought to be. 
These all have the form Xij = xm- The combination of P and the excluded hyperplanes forms 
the vertices of (P, CK), which are all the points of intersection of facets of P and hyperplanes 
in "K that lie in or on the boundary of P. Thus, we count as a vertex every vertex of P itself, 
each point that is the intersection of some facets and some hyperplanes in "K, and any point 
that is the intersection of some hyperplanes and belongs to P, but not intersection points 
that are outside P. (Points of each kind do occur in our examples.) The denominator of 
(P, "K) is the least common denominator of all the coordinates of all the vertices of (P, "K). 
The period of Epo j^(t) divides the denominator; this gives us a known bound on it. 

This geometry might best be explained with an example. Let us consider magic 3x3 
squares, 

Xu Xi2 
X21 X22 X23 
X3I ^32 2:33 

The magic subspace is 

Xu + X12 + Xi3 = X21 + X22 + 3^23 
2 . =2:31+ 2:32 + X33 = Xu + X21 + X31 
= X12+ X22 + X32 = Xis + X23 + X33 
= Xu + X22 + 2:33 = Xi3 + X22 + 2:31 

The hyperplane arrangement [K that captures the distinctness of the entries is 

^ = {(2:11 = 2:12) n s, {xu = X21) n s, . . . , (X32 = X33) n s} . 

Finally, there are two polytopes associated to magic 3x3 squares, depending on whether 
we count them by an upper bound on the entries: 



>o- 



Xu 


X12 


Xl3 




X21 


X22 


X23 


G 


X31 


X32 


X33 





Pe = s n [0, 1]^ 



or by magic sum: 



Pa = s n 



Xu 


X12 


Xl3 




X21 


X22 


X23 


G 


X31 


X32 


X33 





^>0 



2:11 + 2;i2 +Xi3 = l 



Our cubical counting function computes the number of magic squares all of whose entries 
satisfy < Xij < t, in terms of an integral parameter t. These squares are the lattice points 
in 

(p;\U:K)^(iz)^^ 

Our second, affine, counting function computes the number of magic squares with positive 
entries and magic sum t. These squares are the lattice points in 

In general, the number of squares we want to count, N{t), is the Ehrhart quasipolyno- 
mial E°pa r^{t) of an open inside-out polytope (P°,CK). We obtain the necessary Ehrhart 
quasipolynomials by means of the computer program LattE [9j. It computes the closed 
Ehrhart generating function 

Ep(x) := l + ^Ep(t)x* of the values Ep{t) := 4^{P r\ {\zY) . 
t=i 

Counting only interior points gives the open Ehrhart quasipolynomial Epo {t) and its gener- 
ating function 

oo 

Epo(x) := ^Epo{t)xK 

t=i 

Since we want the open inside-out Ehrhart generating function Epo = Ylt^i 

we need several transformations. One is Ehrhart reciprocity [lOj, which is the following 

identity of rational generating functions: 

Epo(t) = (-l)i+^'-^Ep(x-i) . (1) 

The inside-out version [5], Equation (4.6)] is 

E?,o,^(x) = (-l)i+<^'-^Ep,^(x-i). (2) 

We need to express the inside-out generating functions in terms of ordinary Ehrhart gener- 
ating functions. To do that we take the intersection poset 

£(P°, J{) := {P° n nS : S ^ ^} \ {0}, 

which is ordered by reverse inclusion. Note that £(P°,[K) and L{P,'K), defined similarly 
but with P instead of P°, are isomorphic posets because "K is transverse to P; specifically, 
L{P,Ji) = {u : u E £(P°,CK)}, where u is the (topological) closure of u. Now we have the 
Mobius inversion formulas [5l Equations (4.7) and (4.8)] 

Epo,5c(x)= Yl /^(0,«)E„(x) (3) 

and 

Ep,:k(x)= Y1 IMO,m)|E„(x) (4) 
(since is transverse to P; see our general paper [B]). Here fi is the Mobius function of 

£(p°,?c) [n]. 

Thus we begin by getting all the cross-sectional generating functions E„(x) from LattE. 
Then we either sum them by and apply inside-out reciprocity ([2]), or apply ordinary reci- 
procity ([1]) first and then sum by ([3]). (We did whichever of these seemed more convenient.) 
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In the semimagic and magilatin counts we need a third step because the generating func- 
tions we computed pertain to a reduced problem; those of the original problem are obtained 
through multiplication by another generating function. 

Once we have the generating function we extract the quasipolynomial, essentially by the 
binomial series. If an Ehrhart quasipolynomial g of a rational convex polytope has period 
p and degree d, then its generating function q can be written as a rational function of the 
form 

q{x) := q{t) x = -^^ (5) 

7^ (1 - xP) 

for some nonnegative integers oi, 02, . . . , ap(^d+i)- Grouping the terms in the numerator of ([5]) 
according to the residue class of the degree modulo p and expanding the denominator, we 
get 



, , _ Z^r=l 2^j=0 ^PJ+r - _ Y^ 



r=l k>0 



.i=0 



d + k- j 

d 



Hence the rth constituent of the quasipolynomial q is 



Qrit) = apj+r ( P J for r = 1, . . . , p. 



i=o 



d 



2.2. How we apply the method. The initial step is always to reduce the size of the 
problem by applying symmetry. Each problem has a normal form under symmetry, which 
is a strong square. The number of all magic or semimagic squares is a constant multiple of 
the number of symmetry types, because every such square has the same symmetry group. 
For magilatin squares, there are several symmetry types with symmetry groups of different 
sizes, so each type must be counted separately. 

Semimagic and magilatin squares also have an interesting reduced form, in which the values 
are shifted by a constant so that the smallest cell contains 0; and a reduced normal form; the 
latter two are not strong but are aids to computation. Reduced squares are counted either 
by magic sum (the "affine" counting rule) or by the largest cell value (the "cubic" count). 
All reduced normal semimagic squares correspond to the same number of unreduced squares, 
while the different symmetry types of magilatin square give reduced normal squares whose 
corresponding number of unreduced squares depends on the symmetry type. 

The total number of squares, N{t), and the number of reduced squares, R{t), are con- 
nected by a convolution identity N(t) = J2s /(^ ~ where / is a periodic constant (by 
which we mean a quasipolynomial of degree 0; we say constant term for the degree-0 term of 
a quasipolynomial, even though the "constant term" may vary periodically) or a linear poly- 
nomial. Writing for the generating functions N(a;) := Ylt>o ■^i^)-'^'^ ^"^^ similarly R(x), and 
f (^) '■— X]t>o /(^)^*) have N(x) = i{x)Il{x). It follows from the form of the denominator 
in Equation that the period of N divides the product of the periods of / and R. The 
reduced number R{t) is, in the semimagic case, a constant multiple of the number n{t) of 
reduced, normal semimagic squares; in the magilatin case it is a sum of different multiples, 
depending on a symmetry group, of the number of reduced, normal magilatin squares of each 
different type T. Each n{t) is the open Ehrhart quasipolynomial Eqo j{t) of an inside-out 

polytope {Q, J) which is smaller than the original polytope P. 
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We compute n(x) := X]t>o ^('^)'^* from the Ehrhart generating functions Eu(x) of the 
nonvoid sections u of Q° by fiats of J through the following procedure: 

(1) We calculate the fiats and sections by hand. 

(2) We feed each u into the computer program LattE [9J, which returns the closed gener- 
ating function E^(a;), whose constant term equals 1 because u is nonvoid and convex. 

(3) With semimagic squares, by Equations ([l])-(|4l) we have the Mobius-inversion formulas 

n(a;) = E°qo^j{x) = ^/^^(O, u)E„(a;) 

where L := L{Q°,'J), the intersection poset of J). 
The procedure for magilatin squares is similar but taking account of the several types. 

2.3. Checks. We check our results in a variety of ways. 

The degree is the dimension of the polytope, or the number of independent variables in 
the magic-sum equations. 

The leading coefficient is the volume of the polytope. (The volume is normalized so that 
a fundamental domain in the affine space spanned by the polytope has unit volume.) This 
check is also not difficult. The volume is easy to find by hand in the magic examples. In affine 
semimagic the polytope is the Birkhoff polytope B^, whose volume is well known (Section 
14. 2p . The cubical semimagic volume is not well known but it was easy to find (Section I4.1.ip . 
The magilatin polytopes are the same as the semimagic ones. 

The firmest verification is to compare the results of the generating function approach 
with those of direct enumeration. If we count the squares individually for t < ti where 
ti > pd, only the correct quasipolynomial can agree with the counts (given that we know 
the degree d and period p from the geometry). Though ti = pd is too large to reach in 
some of the examples, still we gain considerable confidence if even a smaller value of ti yields 
numbers that agree with those derived from the quasipolynomial or generating function. We 
performed this check in each case. 

2.4. Feasibility. Based on our solutions of the six 3x3 examples we believe our counting 
method is practical. The calculations are simple and readily verified. Linear algebra tells us 
the degree, geometry tells us the period; we obtain the generating function using the Ehrhart 
package LattE \U\ and then apply reciprocity (Equation ([I]) or ([2])) and Mobius inversion 
(Equation (I3l) or (j4j)), and extract the constituents, all with Maple. The programming is not 
too difficult. 

In the magic square problems we found the denominator by calculating the vertices of 
the inside-out polytope. Then we took two different routes. In one we applied LattE and 
Equation ([3]). In the other we calculated N{t) for small values of t by generating all magic 
squares, taking enough values of t that we could fit the quasipolynomial constituents to the 
data. This was easy to program accurately and quick to compute, and it gave the same 
answer. The programs can be found at our "Six Little Squares" Web site [7j. 

In principle the semimagic and magilatin problems can be solved in the same two ways. 
The geometrical method with Mobius inversion gave complete answers in a few minutes 
of computer time after a simple hand analysis of the geometry (see Section H]). Direct 
enumeration on the computer proved unwieldy (at best), especially in the affine case, where 
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the period is largest. A straightforward computer count of semimagic squares by magic 
sum (performed in Maple — admittedly not the language of choice — on a personal computer) 
seemed destined to take a million years. Switching to a count of reduced normal squares, the 
calculation threatened to take only a thousand years. These programs are at our "Six Little 
Squares" Web site [7] , as is a complicated "supernormalized formula" that greatly speeds up 
affine semimagic counting (see Section l4.2.7p . 

2.5. Notation. We use a lot of notation. To keep track of it we try to be reasonably 
systematic. 

M,m refer to magic squares (Section [3]). 
S, s and subscript s refer to semimagic squares (Section Hj). 
L,l and subscript ml refer to magilatin squares (Section |5]). 

R, r refer to reduced squares (the minimum entry is 0), while M, S, L, et al. refer to 
ordinary squares (all positive entries), 
c refers to "cubic" counting, by an upper bound on the entries, 
a refers to "affine" counting, by a specified magic sum. 
X (capital) refers to all squares of that type. 

X (minuscule) refers to symmetry types of squares, or equivalently normalized 
squares. 

3. Magic squares of order 3 
The standard form of a magic square of order 3 is well known; it is 



a + 7 


-a - /3 + 7 


/3 + 7 


-a + /3 + 7 


7 


a - /3 + 7 


-/3 + 7 


a + /3 + 7 


-a + 7 



where the magic sum is s = 37. Taking account of the 8-fold symmetry, under which we may 
assume the largest corner value is a + 7 and the next largest is /3 + 7, and the distinctness 
of the values, we have a > (3 > and a 7^ 2/3. One must also have 7 > a + /3 to ensure 
positivity. 

In this pair of examples, the dimension of the problem is small enough that there is no 
advantage in working with the reduced normal form (where 7 = 0). 

3.1. Magic squares: Cubical count (by upper bound). Here we count by a strict 
upper bound t on the permitted values; since the largest entry is a + /3 + 7, the bound is 
a + /3 + 7 < t. The number of squares with upper bound t is Mc{t). We think of each magic 
square as a t^^-lattice point in P° \ [J'Kc, the (relative) interior of the inside-out polytope 

Pc := {{x,y, z) : < y < X, x + y < z, x + y + z < 1}, 
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[Kc '■= {h} where h : x = 2y, 

but multiplied by t to make the entries integers. Here we use normalized coordinates x = a/t, 
y = f3/t, and z = -f/t. The semilattice of flats is £(P°, :Kc) = {P°,hnP°} with P° < hnP°. 
The vertices are 

O = (0,0,0), C = (0,0,1), D = (l,0,i), i?=(|,|,i), F = (i,l,l), 

of which 0,C,D,F are the vertices of Pc and 0,C,E are those of h (1 Pc. (Both these 
polytopes are simplices.) From Equation ([3]), 

Me(x) = 8E^o,j,^(x) = 8[Epo(x) - Ef,npo{x)] 
which we evaluate by LattE and Ehrhart reciprocity, Equation (JT}: 

- Q R ~ 

^ ^ [(1 -x)2(l-x2)(l-a;4) ' (1 -x)2(l -a;6)_ 

8xi°(2a;2 + 1) 
~ (1 -x)2(l -x4)(l 

_ 8a;^°(2a;2 + l)(a;^ - + l)(a;" + + ■ ■ ■ + a; + l)2(x^° + x^ + ■ ■ ■ + x^ + 1) 

~ (1 • 

From this generating function we extract the quasipolynomial 



Me(t) = 





-16*2+76*- 


-96 _ 


(*-2)(*-6)(*-8) 


if t = 


0,2,6,8 (mod 12) 
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6 ' 




-16*2+73*- 


-58 _ 


(*-l)(*2-15*+58) 


if t = 


1 (mod 12); 




6 




6 ' 




-16*2+73*- 


-102 


(*-3)(*2-13*+34) 


if t = 


3,11 (mod 12); 




6 




6 ' 


< 


-16*2+76*- 


-112 _ 


(*-4)(*2-12t+28) 


if t = 


4,10 (mod 12); 
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6 ' 


t3 


-16*2+73*- 


-90 


(*-2)(*-5)(*-9) 


if t = 


5,9 (mod 12); 




6 




6 ' 


t3 


-16*2+73*- 


-70 _ 


(*-7)(*2-9*+10) 


if t = 


7 (mod 12); 


V 
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6 ' 



and the first few nonzero values for t > 0: 



t 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


M,{t) 


8 


16 


40 


64 


96 


128 


184 


240 


320 


400 


504 


608 


744 


880 


1056 


nicit) 


1 


2 


5 


8 


12 


16 


23 


30 


40 


50 


63 


76 


93 


110 


132 



The last row is the number of symmetry classes, or normal squares, i.e., Mc{t)/8. The rows 
are sequences A108576 and A108577 in the OEIS [I5]. 

The symmetry of the constituents about residue 1 is curious. 

The principal constituent is 

- 16t2 + 76t - 96 _ (t - 2)(t - 6)(t - 8) 
6 ~ 6 ■ 

Its unsigned constant term, 16, is the number of linear orderings of the cells that are induced 
by magic squares. Thus, up to the symmetries of a magic square, there are just two order 
types, even allowing arbitrarily large cell values. (The order types are illustrated in Example 
3.11 of our general magic article, [6].) 
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We confirmed these results by direct enumeration, counting the strongly magic squares 
for t < 60 [7]. 

Compare this quasipolynomial to the weak quasipolynomial: 

(f-l)(t^-2t+3) 



t^-3t^+5f~3 
6 



t3-3t^+8f-6 _ (f-l)(t^-2t+6) 



with generating function 



(x2 + 2x- l){2x 



3 _ ^2 



if t is odd; 
if t is even; 



(1 -X)2(l -X2)2 



3.1.1. Reduced magic squares. A more fundamental count than the number of magic squares 
with an upper bound is the number of reduced squares. Let Rmcit) be the number of 
3x3 reduced magic squares with maximum cell value t, and rmc(^) the number of reduced 
symmetry types, or equivalently of normalized reduced squares with maximum t. Then we 
have the formulas 

t-i t-i 
M,{t) = ^{t-l-k)R,,,,{k) and m,{t) = ^{t - I - k)r^,{k), 

k=0 k=0 

since every reduced square with maximum k gives t — l — k unreduced squares with maximum 
< t (and positive entries) by adding / to each entry where l<l<t — 1 — k. In terms of 
generating functions. 



MJx) 



X 



■ {1-xy 

We deduce the generating functions 



;Rmc{x) and mc(x) 
x^{2x'^ + 1) 



X 



XI 



(9) 



and Rmc(a^ 



'■'"^^"'^ (l-a;4)(l-a;6) 
8rinc(x), and from the latter the quasipolynomial 
' 2t - 16, if t = 
2t-4, 
R^cit) =\ 2t - 8, 
2t 

0, 

1/8-th of these for Vj^dt)) as well as the first few nonzero values: 



if t = 2, 10 
if t = 4, 8 
12, ift = 6 
if t is odd 



> mod 12; 



(10) 



t 


8 


10 


12 


14 


16 


18 


20 


22 


24 


26 


28 


30 


32 


34 


36 


38 


40 


42 


Rmcit) 


8 


16 


8 


24 


24 


24 


32 


40 


32 


48 


48 


48 


56 


64 


56 


72 


72 


72 


^mc(^) 


1 


2 


1 


3 


3 


3 


4 


5 


4 


6 


6 


6 


7 


8 


7 


9 


9 


9 



The sequences are A174256 and A174257 in the OEIS [T5]. 

The principal constituent is 2t — 16, whose constant term in absolute value, 16, is the 
number of linear orderings of the cells that are induced by magic squares — necessarily, the 
same number as with Mc{t). 

We confirmed the constituents by testing them against the coefficients of the generating 
function for several periods. 
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Our way of reasoning, from all squares to reduced squares, is backward; logically, one 
should count reduced squares and then deduce the ordinary magic square numbers from 
them via Equation (Q. Counting magic squares is not hard enough to require that approach, 
but in treating semimagic and magilatin squares we follow the logical progression since then 
reduced squares are much easier to handle. 

3.2. Magic squares: AfRne count (by magic sum). The number of magic squares with 
magic sum t = 37 is Ma(t). We take the normalized coordinates x = a/t and y = (3/t. The 
inside-out polytope is 

Pa: 0<y <x, x + ?/< |, 
!Ka : {h} where h : x = 2y. 

The semilattice of flats is L{P°, JCJ = {P°, h n P°} with P° <hn P°. The vertices are 



O = (0,0), D = a,Q), E 



•2 i\ 

,95 gj? 



6' Q'l 



of which O, D, F are the vertices of Pa and O, E are the vertices of /i fl Pa- From Equations 

©-(ED, 



Ma(x) 



\x 



8[Epo(x) — Efcnp° 



X 



,12 



,12 



;i -X3)2(l 

8x^5(2x3 + 1) 
(l-a;3)(l-x6)(l-s9) 

Sx^^{2x^ + l){x^ + l)(a;i2 



;i -a;3)(l -x9)_ 



x^ + l){x^^ + x^^ ^ Vx^ + 1 



'1 



X 



18^3 



From this generating function we extract the quasipolynomial 



Ma(t) 



2*2 


-32*+ 144 




9 


2*2 


-32t+78 




9 


2*2 


-32t+120 




9 


2*2 


-324+126 




9 


2*2 


-32*+96 




9 


2*2 


-32*+102 


9 



(t2-l6t + 72), ift 



2^+2 
9 



|(t-3)(t-13), 
|(t_6)(t-10). 



(t-7)(t-9). 



i(t-4)(t-12), 
|(t2-16t + 51). 



if t 
if t 
if t 
if t 
if t 



(mod 18); 
3 (mod 18); 
6 (mod 18); 
9 (mod 18); 
12 (mod 18) 
15 (mod 18) 



:iii 



0, 



if t ^ (mod 3); 



and the first few nonzero values for t > 0: 



t 


15 


18 


21 


24 


27 


30 


33 


36 


39 


42 


45 


48 


51 


54 




8 


24 


32 


56 


80 


104 


136 


176 


208 


256 


304 


352 


408 


472 




1 


3 


4 


7 


10 


13 


17 


22 


26 


32 


38 


44 


51 


59 



The last row is the number of symmetry classes, or normalized squares, which is Ma(t)/^ 
The two sequences are A108578 and A108579 in the OEIS tl5j. 
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The principal constituent is 

- 32t + 144 2 



(t^ - 16t + 72), 



9 9 

whose constant term, 16, is the number of hnear orderings of the cells that are induced by 
magic squares — the same number as with Mc{t). 

We verified our results by direct enumeration, counting the strong magic squares for t < 72 

Compare the magic-square quasipolynomial to the weak quasipolynomial: 

'2t2^m+9^ ift = 0(mod3); 
0, ift^0(mod3); 
due to MacMahon [121 Vol. II, par. 409, p. 163], with generating function 

-2x^ + 1 



3^3 



3.2.1. Reduced magic squares. Let Rraa.it) be the number of 3 x 3 reduced magic squares 
with magic sum t, and r^aaii) the number of reduced symmetry types, or equivalently of 
normalized reduced squares with magic sum t. Then 



M.(i) 



i?ma(s) and ma(t) = > 



Q<s<t 
s=t (mod 3) 



0<s<i 
s=t (mods) 



since every reduced square with sum s = t — 3k, where < 3A; < t — 3, gives one unreduced 
square with sum t (and positive entries) by adding 3k to each entry. In terms of generating 
functions, 

„3 



X" 



thus. 



ni^[x) 



[X 



X 



x-^ 



x^\2x'^ + 1) 



(1 -X6)(l -x9) 

and Rnia(a;) = Sr^aix). The quasipolynomial is 



Rma{t) = 8rma(i^) 



3'' 


-16 = 1 


(t-12). 


if t = 


It 


-4 = |( 


t-3). 


if t = 3, 15 


< ¥ 


-8 = |( 


t-6). 


if t = 6, 12 


3'' 


-12 = 1 


(t-9). 


if t = 9 


. 0, 






if t ^ 



> mod 18; 



mod 18. 



(12) 



The initial nonzero values: 



t 


12 


15 


18 


21 


24 


27 


30 


33 


36 


39 


42 


45 


48 


51 


54 


57 


60 


63 


Rmait) 


8 


16 


8 


24 


24 


24 


32 


40 


32 


48 


48 


48 


56 


64 


56 


72 


72 


72 


rmait) 


1 


2 


1 


3 


3 


3 


4 


5 


4 


6 


6 


6 


7 


8 


7 


9 


9 


9 



These sequences are A174256 and A174257 in the OEIS ^ISj. 
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One of the remarkable properties of magic squares of order 3 is that i?inc(2fc) = -Rma(3/c). 
The reason is that the middle term of a reduced 3x3 magic square equals s/3, if s is the 
magic sum, and the largest entry is 2s/3. Thus, the reduced squares of cubic and affine 
type, allowing for the difference in parameters, are the same, and although the counts of 
magic squares by magic sum and by upper bound differ, the only reason is that the reduced 
squares are adjusted differently to get all squares. 

The principal constituent |t — 16 has constant term whose absolute value is the same as 
with all other reduced magic quasipolynomials. 

We confirmed the constituents by comparing their values to the coefficients of the gener- 
ating function for several periods. 

4. Semimagic squares of order 3 

Now we apply our approach to counting semimagic squares. Here is the general form of a 
reduced, normalized 3x3 semimagic square, in which the magic sum is s = 2a + 2/3 + 7: 






/3 


2a + /3 + 7 


a + /3 


a + (3 + -f-6 


6 


a + /3 + 7 


a + 6 


13-5 



(13) 



Proposition 1. A reduced and normalized 3x3 semimagic square has the form f[T3|) with 
the restrictions 

< a, /3, 7; 

(14) 



and 



/3+a+7. 




(15) 



The largest entry in the square is w := X13 = 2a + /3 + 7. 

Each reduced normal square with largest entry w corresponds to exactly 72(t — w — 1) 
different magic squares with entries in the range (0,t), for < w < t. Each reduced normal 
square with magic sum s corresponds to exactly 72 different magic squares with magic sum 
equal to t, if t = s (mod 3), and none otherwise, for < s < t. 

Proof. By permuting rows and columns we can arrange that Xn = minxij and that the top 
row and left column are increasing. By flipping the square over the main diagonal we can 
further force X21 > Xu- By subtracting the least entry from every entry we ensure that 
Xii = 0. Thus we account for the 72{t — w — 1) semimagic squares that correspond to each 
reduced normal square. 
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The form of the top and left sides in f|T3l) is explained by the fact that xu < Xi2 < 
^21 < < X13. The conditions Xij > Xu for i,j = 2,3, together with the row-sum and 
column-sum equations, imply that X13 is the largest entry and that X23 < Xi2- 

The only possible equalities amongst the entries are ruled out by the following inequations: 

3^22 7^ 3^12, 3^21, 2:23; 
3^32 7^ X12, X22] 
2^33 7^ 3:23) 3:32. 

These correspond to the restrictions f|T5|) : 

2:22 7^ X12 
X22 7^ X21 

X22 7^ a;23 
X32 ^ X12 

X32 7^ X22 
7^ a;23 
2^33 7^ 3:32 

4.1. Semimagic squares: Cubical count (by upper bound). We are counting squares 
by a strict upper bound on the allowed value of an entry; this bound is the parameter t. Let 
Sc{t), for t > 0, be the number of semimagic squares of order 3 in which every entry belongs 
to the range (0, t). 

4.1.1. Counting the weak squares. The polytope Pc is the 5-dimensional intersection of [0, 1]^ 
with the semimagic subspace in which all row and column sums are equal. This polytope is 
integral because it is the intersection with an integral polytope of a subspace whose constraint 
matrix is totally unimodular; so the quasipolynomial is a polynomial. By contrast, the 
inside-out polytope (P, !K) for enumerating strong semimagic squares has denominator 60. 
We verified by computer counts for t < 18 that the weak polynomial is 

3t^ - ISt^ + 35^3 - 45^2 + 32t - 10 _ (t - l){t^ - 2t + 2){3t'^ - 6t + 5) 
10 ~ 10 

with generating function 

{7x^ -2x + l){2x^ + + 4x - 1) 
(1 ■ 

We conclude that Pc has volume 3/10. 

4.1.2. Reduction of the number of strong squares. We compute Sc{t) via Rc{w), the number 
of reduced squares with largest entry equal to w. The formula is 

t-i 

Sc{t) = J2(t~'^-w)Rc{w). (16) 

w=0 
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5 7^ a + 7; 

^+« + 7 
^ 2 

6 (3 — a; 



/3 + 7 



6^a + 6y^ 13-6 



13 — a 



The value Rc{w) = 72rc{w) where rc{w) is the number of normal reduced squares (in which 
we know the largest entry to be X13). Thus r^ls) counts the number of ^-integral points in 
the interior of the 3-dimensional polytope Qc defined by 

0<x,y; 0<z<y; x + y<^ (17) 



with the seven excluded (hyper)planes 

y — X y 1 — y — 2x 1 — x — y 



y — X, 1 — a; — 2y, 1 — 2x — 2y 



2 ' 2' 2 ' 2 

the three coordinates being x = a/w, y = (3/w, and z = 6/w. 

The hyperplane arrangement for reduced normal squares is that of (|T8i) . We call it Jc- 
Thus re(s) = E^o^j^{s). 

4.1.3. Geometrical analysis of the reduced normal polytope. We apply Mobius inversion, 
Equation (E]), over the intersection poset L{Q°,Jc)- We need to know not only L{Q°,'Jc) 
but also all the vertices of (Qc'Jc), since they are required for computing the Ehrhart gen- 
erating function and estimating the period of the reduced normal quasipolynomial. 
We number the planes: 





X — y + 2z = 


0, 




y-2z = 


0, 




2x + 2z = 


1, 


7r4 


x + 2z = 


1, 




X — y + z = 


0, 




X + y + z = 


1, 




2x + y + z = 


1. 



The intersection of two planes, ttj fl vTfc, is a line we call Ijk, vts fl tts fl ttq is a line we also 
call The intersection of three planes is, in general, a point but not usually a vertex of 

Our notation for the line segment with endpoints X, Y is XY, while XY denotes the 
entire line spanned by the points. The triangular convex hull of three noncollinear points 
X, Y, Z is XYZ. We do not need quadrilaterals, as the intersection of each plane with Q is 
a triangle. 

We need to find the intersections of the planes with Q°, separately and in combination. 
Here is a list of significant points; we shall see it is the list of vertices of {Qc, "Jc)- The first 
column has the vertices of Qc, the second the vertices of (Qc, Jc) that lie in open edges, the 
third the vertices that lie in open facets, and the last is the sole interior vertex. 



(19) 



The denominator of {Qc, Jc) is the least common denominator of all the points; it evidently 
equals 4 ■ 3 ■ 5 = 60. 
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= 


(0,0,0), 


Dc = 


(0,i,i)GOC, 


Fc = 


(0, 


2 

3' 


\) 


G DBG, He = (|, |, 1 


A = 


(io,o). 


Ec = 


{\,\,0)eAB, 


Gc = 


(i 


1 

2' 


\) 


G ABG, 


Bc = 


(0,1,0), 


K = 


{\, \, \) e AC, 


G'c = 




3 

5' 


\) 


G ABG, 


Cc = 


(0,1,1), 


e: = 


(0,l,|) e5C, 


g: = 


(i 


3 

5' 


f) 


G ABG, 



Plane 




Intersection with edge 




Intersection 




OA 


OB 


OC 


AB 


AC 


BC 


with Q 













E 


iQ 


E" 


OEE" 




OA 








A 


A 


E" 


CAE" 




A 





D 


A 


A 


E" 


ADE" 


7r4 


iQ 





D 


iQ 


E' 


E" 


DE"E' 










OC 


E 


C 


C 


OCE 






B 


D 


B 


E' 


B 


BDE' 




A 


B 


D 


AB 


A 


B 


ABD 



Table 1. Intersections of planes of J with edges of Q. A vertex of Q contained 
in TTj will show up three times in the row of VTj. In order to clarify the geometry, 
we distinguish between a plane's meeting an edge line outside Q and not 
meeting it at all (i.e., their being parallel). 



The intersections of the planes with the edges of Qc are in Table [T] The subscript c is 
omitted. 

Table |2] shows the lines generated by pairwise intersection of planes. 





1^2 




Ha 




7r6 






s = 
y = 2z 


X + Z = \ 

y-z = \ 


y = l 
x + 2z = 1 


z = 
X = y 


z = 2y-l 
X = 2 — 3y 


X + z = ^ 

y-z = \ 


7r2 




y = 2z 
X + z = \ 


y = 2z 
X + y = 1 


y = 2z 

X = z 


y = 2z 
X + 3z = 1 


y = 2z 
2x + 3z = 1 








X = 
z = ^ 

^ 2 


X + z = \ 
y = 2 


y = z 










X = l-2z 
y = 1- z 


x = l-2z 
y = z 


X = l-2z 
y = 3z-l 












^356 


X = 1 — 2y 
z = 3y-l 














X = 
y + z = 1 



Table 2. The equations of the pairwise intersections of planes of Jc- 



In Table [3] we describe the intersection of each line with Qc and with its interior. The 
subscript c is omitted. 

Last, we need the intersection points of three planes of Jc', or, of a plane and a line. Some 
are not in Qc at all; them we can ignore. Some are on the boundary of Qc', they are necessary 
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vr2 










Try 






E" 
(BC) 


E" 
(BC) 


OE 
(OAB) 




EE 






AE" 
{ABC) 


E" 

(BC) 


OG 


EG 


AE 








DE" 
(OBC) 


^356- DG 


AD 
{OAG) 


7r4 








DG" 


DE' 
(OAG) 


D 

(OG) 












^356 


ED 














BD 
{DBG) 



Table 3. The intersections of lines with Q and Q° . The second (parenthe- 
sized) row in each box shows the smallest face of Q to which the intersection 
belongs, if that is not Q itself; these intersections are not part of the intersec- 
tion poset of (Q°, J). 



in finding the denominator, but all of them are points already listed in f lT9|) . It turns out 
that 



7r2 n TTs n TTy = He 



is the only vertex in Q°^, so it is the only one we need for the intersection poset. 

Here, then, is the intersection poset (Figured]). The subscript c is omitted. In the figure, 
for simplicity, we write vTj, etc., when the actual element is the simplex vTj fl Q° , etc.; we also 
state the vertices of the simplex. The Mobius function /i(0, u) equals (_i)c°dim« ^i^e 
exception of /i(0, I^^q) = 2. 



4.1.4. Generating functions and the quasipolynomial. That was the first half of the work. 
The second half begins with finding rc{w) = E°qo j -)(u') from the Ehrhart generating func- 
tions Eu{t) of the intersections by means of ([6]). The next step, then, is to calculate all 
necessary generating functions. This is done by LattE. r^ix) is the result of applying reci- 
procity to the sum of all these rational functions (with the appropriate Mobius-function 
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H 



k5 {DG") he {FG') Zi7 {EF) k^e (DG) (FG) h, {ED) I25 (OG) h, (AF) 




7r4 {DE"E') TTi {OEE") ttb {BDE') tts {OCE) tts {ADE") ttt {ABD) tts {OAE") 





W (OABC) 

Figure 1. The intersection poset L{Q°,3) for semimagic squares. The dia- 
gram shows both the fiats and (in parentheses) their intersections with Q. 



multipher —1, excepting Edg°{x) whose multiplier is —2). The result is: 

-rc(l/x) = BoABcix) +'EoEE"{x) +'EoAE"{x) 

+ 'E'ADE"{x) + 'Ede'e"{x) + 'Eoce{x) 

+ 'EbDE'{x) + BABoix) + BFG'ix) 

+ Eef{x) + Eog{x) + Efg{x) 

+ Baf{x) + 2I^dg(x) + ^DG"(x) 

+ Bde(x) + Bh(x) 

1 1 



+ 



+ 



{l-xf{l-x^) {1 - x){l - x^){l - x^) {l-x){l-x^y 

1 1 1 

+ (1 _ ^2)3 + (1 _ ^2)2(1 _ ^3) + (1 _ ^yi^l _ ^3-) 

1 11 

I \ / -I o\/-. o\ /-I \ / -I o\o 



{1 - x){l - x^){l - x^) {l-x){l-xY (l-a;3)(l-a;5) 



+ 



+ 



+ 



(1-X3)2 (l-x)(l-X^) (1-X^)(1-X^) 

111 

+ VZ KTVZ ^ + 2- — + 



+ 



(1-X2)(1-X3) (l-x2)(l-x4) {1-X^){1-X^) 

1 1 



{1 — x'^){l — x^) 1 — x^ 

18 



(20) 



Then by f fT6|) the generating function for cubically counted semimagic squares is 
Sc(x) = 72- r^rc(x) 

_ 72x1° [18 + 46 + 69 x^ + 74 x^ + 65 + 46 + 26 a;^ + 11 + 4 x + 1] 

~ (l-x2)2(l-a;3)2(l-a;4)(l-x5) ' 

(21) 

From the geometrical denominator 60 or the (standard-form) algebraic denominator (1 — 
^,60 -J 5 know the period of Sdt) divides 60. We compute the constituents of Sdt) by the 
method of Section I2.lt the resuh is that 



10 
3 



75 



75 



331 
331 



UO 8 3 
where Ci varies with period 6, given by 

1464, 



5989 
~W 

11933 
20 



+ Ci{t)t — co(t), if t is even; 



+ Ci(t)t — Co(t), if t is odd; 



(22) 



Ci(t) 



(mod 6); 



if t = 0,2 
if t = 4 
if t = 1 
if t = 3, 5 

and Co, given by Table HJ varies with period 60. (It is curious that the even constant terms 
have half the period of the odd terms.) Thus the period of Sc turns out to be 60, the largest 
it could be. 

The principal constituent (for t = 0) is 



1456, 

2831 

2 ' 
2847 

2 ' 



3 = 75 4 331 o 

TTT* - -irt + — t 



5989 2 , ... 
-t^ + 1464t 



1296. 



10 8 3 10 

Its unsigned constant term, 1296, is the number of order types of semimagic squares. Al- 
lowing for the 72 symmetries of a semimagic square, there are just 18 symmetry classes of 
order types. 

For the first few nonzero values of Sc{t) see the following table. (This is sequence Al 73546 
in the OEIS [IS].) The third row is the number of normalized squares, or symmetry classes 
(sequence A173723), which equals Sdt) /72. The other lines give the numbers of reduced 
squares (sequence A173727) and of reduced normal squares (i.e., symmetry types of reduced 
squares; sequence A173724), which may be of interest. 



t 


8 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


Sc{t) 








72 


288 


936 


2592 


5760 


11520 


20952 


35712 


57168 


88272 


Sc{t) 








1 


4 


13 


36 


80 


160 


291 


496 


794 


1226 


Rc{t) 


72 


144 


432 


1008 


1512 


2592 


3672 


5328 


6696 


9648 


11736 


15552 


rdt) 


1 


2 


6 


14 


21 


36 


51 


74 


93 


134 


163 


216 



Compare the strong to the weak quasipolynomial. The leading coefficients agree and 
the strong coefficient of is constant. These facts, of which we made no use in deducing 
the quasipolynomial, provide additional verification of the correctness of the counts and 
constituents. 
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t 


co(t) 


t 


co(t) 


t 


co(t) 


t 


co(t) 


t 


co(t) 





1296 


12 


1296 


24 


6624 
5 


36 


6192 
5 


48 


6624 
5 


1 


110413 
120 


13 


120781 
120 


25 


23465 
24 


37 


23465 
24 


49 


120781 
120 


2 


3824 
3 


14 


19552 
15 


26 


18256 
15 


38 


19552 
15 


50 


3824 
3 


3 


47727 
40 


15 


9315 
8 


27 


9315 
8 


39 


47727 
40 


51 


44271 
40 


4 


18152 
15 


16 


16856 
15 


28 


18152 
15 


40 


3544 
3 


52 


3544 
3 


5 


25705 
24 


17 


25705 
24 


29 


131981 
120 


41 


121613 
120 


53 


131981 
120 


6 


6192 
5 


18 


6624 
5 


30 


1296 


42 


1296 


54 


6624 
5 


7 


25193 
24 


19 


129421 
120 


31 


119053 
120 


43 


129421 
120 


55 


25193 
24 


8 


19552 
15 


20 


3824 
3 


32 


3824 
3 


44 


19552 
15 


56 


18256 
15 


9 


44847 
40 


21 


41391 
40 


33 


44847 
40 


45 


8739 
8 


57 


8739 
8 


10 


3544 
3 


22 


3544 
3 


34 


18152 
15 


46 


16856 
15 


58 


18152 
15 


11 


130253 
120 


23 


140621 
120 


35 


27433 
24 


47 


27433 
24 


59 


140621 
120 



Table 4. Constant terms (without the negative sign) of the constituents of 
the semimagic cubical quasipolynomial Sdt). 



4.1.5. Another method: Direct counting. We checked the constituents by directly counting 
(in Maple) all semimagic squares for t < 100. The numbers agreed with those derived from 
the generating function and quasipolynomial above. 

4.2. Semimagic squares: AfRne count (by magic sum). Now we count squares by 
magic sum: we compute 5'a(t), the number of squares with magic sum t. 

4.2.1. The Birkhojf polytope. The polytope P for semimagic squares of order 3, counted by 
magic sum, is 4-dimensional and integral. (It is the polytope of doubly stochastic matrices 
of order 3, i.e., a Birkhoff polytope [Hill].) 

4.2.2. Affine weak semimagic. The polytope for weak semimagic squares of order 3 is the 
same P. 

The weak quasipolynomial, or rather, polynomial, first computed by MacMahon [T2^, Vol. 
II, par. 407, p. 161], is 

t^ - 6t^ + 15t2 _ I8t + 8 _ (t - l)(t - 2)(t2 - 3t + 4) 
8 ~ 8 

with generating function 

Qx^ - 9x^ + lOx^ - 5x + 1 
(1 -x)5 ■ 
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4.2.3. Reduction. The count is via Rs,{s), the number of reduced squares with magic sum s. 
The formula is 



if t > 0. 



(23) 



0<s<t-3 
s=t (mods) 



We have Ra{s) = 72ra(w), where ra(s) is the number of reduced, normahzed squares with 
magic sum s, equivalently the number of ^-integral points in the interior of the 3-dimensional 
polytope Qa defined by 

0<x,y; 0<z<y] x + y<l- (24) 



with the seven excluded (hyper)planes 

y — X y 1 — y — 2x 1 — x — y 



, y — X, 1 — X — 2y, 1 — 2x — 2y 



(25) 



2 ' 2' 2 ' 2 
the three coordinates being x = a/s, y = (3/ s, and z = 6/s. 

The hyperplane arrangement for reduced, normalized squares is that of fl25|l . We call it 
Ja- Thus ra(s) = E°o jjs). 

4.2.4. The reduced, normalized weak polytopal quasipolynomial. This function simply counts 
i-lattice points in Q^. The counting formula is X^qX^/sS^I; summed over all triples that 
satisfy f|T4|) . It simplifies to 



which gives the Ehrhart quasipolynomial 

s— 
2 

3 



2 



s-2 
2 

3 



^(s — l)(s — 3)(s — 5), for odd s; 
^(s — 2)(s — 4)(s — 6), for even s. 



(26) 



The leading coefficient is vol 
We deduce from fl26l) that 



and by reciprocity that 



Epo(x) 



Ep(x) 



x''(l + x) 
(l-a;2)4 

1+x 



'1 — X 



2^4' 



4.2.5. Geometrical analysis of the reduced, normalized polytope. We apply Mobius inversion. 
Equation ([6]), over the intersection poset L{Ql,J). 
We number the planes: 



7r4 



X — ?/ + 2^; = 0, 
y-2z = 0, 
2x + y + 2z = 1, 
X + y + 2z = 1, 
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TTs : X — y + z = 0, 
ttq : X + 2y + z = 1, 
ir-j : 2x + 2y + z = l. 

The intersection of two planes, iTj fl tt^, is a line we call Ijk', tts fl tts fl ttq is a line we also 
call /sse- The intersection of three planes is, in general, a point but not usually a vertex of 
(Qa, Ja)- Our geometrical notation is as in the cubical analysis. 

We need to find the intersections of the planes with Q°, separately and in combination. 
Here is a list of significant points; we shall see it is the list of vertices of (Qa, Ja)- The first 
column has the vertices of Qa, the second the vertices of (Qa, ^a) that lie in open edges, the 
third the vertices that lie in open facets, and the last is the sole interior vertex. 



= 


(0,0,0), 






\) e oc, 


Fa = 






G OBC, Hs. = (y, |, y 


A = 


(io,o). 




V4' 4' 


0) G AB, 


Ga = 


V6' 6' 


1) 


G ABC, 




(0,|,0), 


EL = 


U' 4' 


\) e AC, 


G'a. = 


{- ^ 

V8' 8' 


1) 


G 


Ca = 


fo i 


K = 


(o,i 


\) e BC, 


f-^a — 


f- ^ 

V8' 8' 


1) 


G AEC. 



The least common denominator of 0,A,B^,C^ explains the period 2 of -Eq". The denom- 
inator of (Qa, Ja) is the least common denominator of all the points; it evidently equals 
8 . 3 ■ 5 ■ 7 = 840. 

The intersections of the planes with the edges of Q^, are in Table [T] Table E] shows the 
lines generated by pairwise intersection of planes. Table [3] describes the intersection of each 
line with and with Q°^. 





7r2 


7r3 








7l7 




X = 
y = 22; 


X - 

y = '^ 


x + 2z = \ 

y = \ 


X = y 
z = 


X = 2 — 5y 
z = 3y-l 


rr — 1-52; 
X — 4 

y - 






x + y = \ 
y = 2z 


X = 1 — 2y 
y = 2z 


X = z 
y = 2z 


X = 1 — 5z 
y = 2z 


rr — l-5z 
X — 2 

y = 2z 








X = 

y = l-2z 


X + z = ^ 
'356 : 1 
^=3 


2x + 3y = 1 
z = y 










X = 3y — 1 
z = l-2y 


X = 1 — 3y 
z = y 


x + y = I 
z = ^ 

^ 3 












^356 


X = 1 — 3y 
z =Ay-l 














X = 

z = l-2y 



Table 5. The equations of the pairwise intersections of planes of J, 



Last, we need the intersection points of three planes of Ja; or, of a plane and a line. Some 
are not in Qa at all; them we can ignore. Some are on the boundary of Qa; they are necessary 
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in finding the denominator, but all of them are points already listed in (1271) . It turns out 
that 

7r2 n 775 n 7r7 = 

is the only vertex in Qa, so it is the only one we need for the intersection poset. 

The combinatorial structure and the intersection poset (Figure [1]) for the affine count are 
identical to those for the cubical count. The reason is that the affine polytope Pa is the 
4-dimensional section of by the flat in which the magic sum equals 1, and this fiat is 
orthogonal to the line of intersection of the whole arrangement 'K^.. 



4.2.6. Generating functions and the quasipolynomial. The second half of the affine solution is 
to find ra(s) = E°qo 3^){^) by applying Equations after finding the Ehrhart generating 

functions E„(s) for u G L{Ql,3a)- The next step, then, is to calculate those generating 
functions. This is done by LattE. Then (— l)^ra(x~^) is the sum of all these rational functions; 
that is, 

-ra(l/x) = EoABc{x) + EoEE"{x) + EoAE"{x) 

+ Eade"{x) + Ede"E'{x) + Eoce{x) 

+ EBDE'ix) + EABoix) + EpG'ix) 

+ Eef{x) + Eog{x) + Efg{x) 
+ Eaf{x) + 2Edg{x) + EDG"ix) 
+ Ede{x) + Eh{x) 

1 1 1 



(l-x)(l-x2)3 (1 _a;)(l (1 -x)(l -a;2)(l -x^) 

1 1 1 

+ 
+ 



(1 -x2)(l -a;3)(l-x4) (1 _a;3)(l -x4)2 (l - x){l - x^){l - x^) 

1 1 1 



{1 - x^){l - x^){l - x^) (1 - a;2)2(l - a;3) {1 - x^){l - x^ 
1 1 1 

+ TZ TTTZ ^ + 



(l-x4)(l-x5) (1-X)(1-X6) (l-x5)(l-x6) 

1 1 1 
+ 7Z ^TT: ^ + 2- ^-7- ^ + 



+ 



(1 -x2)(l (i_a;3)(l-x6) {1 - x^){l - x^) 
1 1 



(1 -x3)(l -x^) 1 -x^ ■ 
The generating function for the affine count of semimagic squares, by ( l23l) . is 



(28) 



x^ 



Sa(a;) = 72 ra(a;) 

1 — x-^ 



72x 



15 



18x^^ + 5x^° + 15x^^ + llx^^ - 8x^^ + x^^ - 23x^^ - 13x^^ - 22x^^ - 9x"" 
- 16x^° + - 3x® + 7x^ + 7x^ + 9x^ + 7x^ + 6x^ + 4x2 + 2x + 1 

(1 - X3)2(l - x4)(l - x5)(l - x6)(l - x7)(l - x8) 
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(29) 



From the geometrical or generating-function denominator we know that the period of 
S'a(t) divides 840 = lcm(3, 4, 6, 7, 8). This is long, but it can be simplified. The factor 7 in 
the period is due to a single term in ( l28l) . If we treat it separately we have ra as a sum 
of the H-teim /{I — x^) and a "truncated" generating function for v^{x) + jz^^ and a 
corresponding truncated expression 

72- 



Sa(a;) 



.7^ 



-72x 



10 



llx^^ + 5x^^ + 15x^^ + x"° + I2x 



„18 

,11 



,16 



15 



iX 



9x^° - - - - + + - 1 



(1 - X3)2(l - - x5)(l - x6)(l - X^) 



We extract the constituents from this expression as in Section I2.H separately for the two 
parts of the generating function. The constituents are all of the form 



- + a^{t)t^ - ai{t)t + ao{t) - 72Sr{t), 
8 2 



where 5*7 (t) is a correction, to be defined in a moment, and 



r^, ift 



^, ift 
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1,5 
2,4 
if t = 3 



(mod 6); 



(30) 



ai{t) 



( 1968 
5 

1158 
5 

1383 
5 

1653 
5 

1428 
5 

1923 
5 

1113 
5 

1698 
5 



if t = 

if t = 1,5 
if t = 2, 10 
if t = 3 
if t = 4, 8 
if t = 6 
if t = 7, 11 
if t = 9 



(mod 12); 



and ao(t) is given in Table El 

We call the constituents of the quasipolynomial 

5a(t) + 12S^{t) = - + a2{t)e - ai{t)t + ao(t) 
o 2 

the truncated constituents of Sa{t), since they correspond to the truncated generating func- 
tion mentioned just above. The term that undoes the truncation is 



t-1 
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1, if t= 10,13, 16,17, 19,20 (mod 21); 
0, otherwise 



t-t 
21 



24 



where t :— the least positive residue of t modulo 7 and 

^ _ (l, ift= 10, 13, 16, 17, 19, 20 (mod 21); 
1 0, otherwise. 

Note that f = 21 if i = 0, so that ^7(0) = -1 and in general S7{21k) = k - 1. 



t 


ao{t) 


t 


ao{t) 


t 


ao(t) 


t 


ao(t) 


t 


ao{t) 


t 


ao(t) 





1224 


20 


524 


40 


584 


60 


1188 


80 


560 


100 


548 


1 


7259 


21 


31419 
/in 


41 


6299 

/in 
4U 


61 


8699 

/in 
4U 


81 


29979 

/in 
4U 


101 


7739 

/in 
4U 


2 


1801 
5 


22 


1741 

5 


42 


5121 

5 


62 


1621 

5 


82 


1921 
5 


102 


4941 

5 


3 


23067 
40 


23 


827 
40 


43 


347 
40 


63 


24507 
40 


83 


-613 
40 


103 


1787 
40 


4 


2452 




24 


5832 




44 


2332 




64 


2632 




84 


5652 




104 


2512 




5 


2239 

Q 


25 


2143 

Q 


45 


6975 

Q 


65 


1951 

Q 


85 


2431 

Q 


105 


6687 

Q 


6 


4653 




26 


1513 




46 


1453 




66 


4833 




86 


1333 



106 


1633 




7 


5243 

/I A 


27 


26523 

/in 
4U 


47 


4283 

/in 
4U 


67 


3803 

/in 
4U 


87 


27963 

/in 
4U 


107 


2843 

A n 
4U 


8 


2224 
5 


28 


2164 
5 


48 


5544 
5 


68 


2044 
5 


88 


2344 
5 


108 


5364 
5 


9 


31131 
40 


29 


8891 
40 


49 


8411 
40 


69 


32571 
40 


89 


7451 
40 


109 


9851 
40 


10 


413 


30 


1017 


50 


389 


70 


377 


90 


1053 


110 


353 


11 


539 

A n 
4U 


31 


2939 

/in 
4U 


51 


24219 

/in 
4U 


71 


1979 

/in 
4U 


91 


1499 

/in 
4U 


111 


25659 

A n 
4U 


12 


5796 




32 


2656 




52 


2596 




72 


5976 




92 


2476 




112 


2776 

re 




13 


7547 

40 


33 


28827 
40 


53 


6587 
40 


73 


6107 
40 


93 


30267 
40 


113 


5147 
40 


14 


1477 
5 


34 


1777 
5 


54 


4797 
5 


74 


1657 
5 


94 


1597 
5 


114 


4977 
5 


15 


5823 
8 


35 


799 
8 


55 


1279 
8 


75 


5535 
8 


95 


1087 
8 


115 


991 

8 


16 


2488 
5 


36 


5508 
5 


56 


2368 
5 


76 


2308 
5 


96 


5688 
5 


116 


2188 
5 


17 


8603 
10 


37 


11003 


57 


32283 


'I'i 


10013 


97 


9563 
10 


117 


33723 


10 


10 


10 


40 


18 


4689 
5 


38 


1189 
5 


58 


1489 
5 


78 


4509 
5 


98 


1369 
5 


118 


1309 
5 


19 


2651 
40 


39 


26811 
40 


59 


1691 
40 


79 


4091 
40 


99 


25371 
40 


119 


3131 
40 



Table 6. Constant terms of the truncated constituents of Ss,{t). 



The period of the constant term of the truncated constituents is 120. It follows that 
has period 7 • 840, that is, 840. 

The principal constituent of Sa,{t) (that is, for t = 0) is 



9 . 



243 



25 



13896 
"35" 



-t + 1296. 



(This incorporates the effect of the term —725*7.) The constant term is the same as in the 
cubic count, as it is the number of order types of semimagic squares. 

We give the first few nonzero values of S'a(t) in the following table. (This sequence is 
Al 73547 in the OEIS [ISj.) The third row is the number of normalized squares, or symmetry 
classes (sequence A173725); this is S'a(t)/72. The last rows are the numbers of reduced 
squares (sequence Al 73728) and of reduced, normahzed squares (sequence A173726) with 
magic sum t. 



t 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


5a(t) 











72 


144 


288 


576 


864 


1440 


2088 


3024 


3888 


5904 


Sa(t) 











1 


2 


4 


8 


12 


20 


29 


42 


54 


82 




72 


144 


288 


504 


720 


1152 


1512 


2160 


2448 


3816 


3960 


5544 


6264 




1 


2 


4 


7 


10 


16 


21 


30 


34 


53 


55 


77 


87 



4.2.7. Alternative methods: Direct counting and direct computation. We verified our formulas 
by computing Sg,{t) for t < 100 through direct enumeration of normal squares. The results 
agree with those computed by expanding the generating function. 

We also applied Proposition [1] to derive a formula, independent of all other methods, by 
which we calculated numbers (which we are not describing; see the "Six Little Squares" Web 
page ^) that allowed us to find the 840 constituents by interpolation. These interpolated 
constituents fully agreed with the ones given above. 

5. Magilatin squares of order 3 

A magilatin square is like a semimagic square except that entries may be equal if they 
are in different rows and columns. The inside-out polytope is the same as with semimagic 
squares except that we omit those hyperplanes that prevent equality of entries in different 
rows and columns. Thus, in our count of reduced squares, we have to count the fractional 
lattice points in some of the faces of the polytope. 

The reduced normal form of a magilatin square is the same as that of a semimagic square 
except that the restrictions are weaker. It might be thought that this would introduce 
ambiguity into the standard form because the minimum can occur in several cells, but it 
turns out that it does not. 

Proposition 2. A reduced, normal 3x3 magilatin square has the form f[T^ with the re- 
strictions 

0</3,7; 

< «; (31) 
< 5 < /3; 

and (fT5|) . Each reduced square with w in the upper right corner corresponds to exactly 
t — w — 1 different magilatin squares with entries in the range {0,t), for < w < t. Each 
reduced square with magic sum s corresponds to one magilatin square with magic sum equal 
to t, if t = s (mod 3), and none otherwise, for < s < t. 

Proof. The proof is similar to that for semimagic squares; we can arrange the square by 
permuting rows and columns and by reflection in the main diagonal so that Xu is the 
smallest entry, the first row and column are each increasing, and X21 > 2:12. We cannot say 
X21 > X12 because entries that do not share a row or column may be equal. Still, we obtain 
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the form f|T3l) with the bounds fl3T|) and the same inequations f|T5|) as in semimagic because 
all the latter depend on having no two equal values in the same line (row or column). ■ 



Each reduced, normal magilatin square gives rise to a family of true magilatin squares by 
adding a positive constant to each entry and by symmetries, which are generated by row 
and column permutations and reflection in the main diagonal. Call the set of symmetries 
&. As with semimagic, |C5| = 2(3!)^ = 72. Each normal, reduced square S gives rise to 
l^/^sl = 72/\&s\ squares via symmetries, where ©5 is the stabilizer subgroup of S. If all 
entries are distinct, then the square is semimagic, (3s is trivial, and everything is as with 
semimagic squares. However, ifQ; = Oor(5 = Oor5 = /3, the stabilizer is nontrivial. We 
consider each case in turn. 

The case a = < 6 . Here 6 < P because no line can repeat a value. To fix the square we 
cannot permute any rows or columns but we can reflect in the main diagonal, so \&s\ =2. 
Moreover, (113]) reduces to 

We are in OBC, the x = facet of Q, with the induced arrangement of three lines, J^^°. 
The number of reduced magilatin squares of this kind is 36roBcit), where roBcif) is the 
number of ^-lattice points in the open facet and, equivalently, the number of symmetry types 
of reduced magilatin squares of this kind. We apply Equations ([I])-(l4]) to the intersection 
poset £(0-BC°, which is found in Figure El The Mobius function fi{OBC,u) equals 

^ codim u 

The case 6 = < a. In this case a nontrivial member of ^5 can only exchange the two 
zero positions. Such a symmetry that preserves the increase of the first row and column is 
unique (as one can easily see); thus = 2. Furthermore, (fT5|) reduces to 

We are in OAB, the z = facet, with the induced arrangement J^^" of one line. The number 
of reduced magilatin squares of this kind is SGroAsit), where roAsit) is the number of ~ 
lattice points in the open facet and, equivalently, the number of symmetry types of reduced 
squares. We apply Equations to the intersection poset L{0AB°,3^^^), shown in 

Figured The Mobius function fi{OAB,u) equals (-i)™dim« 

The case 6 = f3. Here we must have a > 0. There are two zero positions in opposite 
corners. A symmetry that exchanges them and preserves increase in the first row and column 
is uniquely determined, so \(3s\ = 2. The inequations reduce to 

/3 7^ 7,« + 7- 

We are in OAC, the facet where y = z, with the two-line induced arrangement J^"^. The 
number of reduced magilatin squares of this kind is 36roAc{t), where roAcit) is the number 
of i-lattice points in the open facet, equally the number of reduced symmetry types. We 
apply Equations to the intersection poset L{OAC° ,3^^'') in Figured The Mobius 

function fi{OAC,u) equals (-i)™dim«_ 

The case a = = S. In these squares there are three zero positions and the whole square 
is a cyclic latin square. Any symmetry that fixes the zero positions also fixes the rest of 
the square. There are 3! symmetries that permute the zero positions, generated by row 
and column permutations. They all preserve the entire square. Therefore \&s\ = 6. The 
inequations disappear. We are in the edge OB, which is the face where x = z = 0, with the 
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OE 



OAB 




OBC 

L {OAB\ a^=0) L {OAC\ Oy=') L {OBC°, a^=0) 



G' G G' 




£(^5C°, EE" BE' AE" GE E'E 




ABG 

Figure 2. The four facet intersection posets of {Q°,J). 



empty arrangement, J^=^=o = 0. The number of reduced magilatin squares of this kind is 
12roB(t), where rosif) is the number of —lattice points in the open edge, also the number 
of reduced symmetry types. The intersection poset L{OB° , 0) consists of the one element 
OB, whose Mobius function n{OB,OB) = 1. 

To get the intersection posets we may examine Tables [U and |3] to find the edges and vertices 
of {Q°,J) in each closed facet. We also need to know which vertex is in which edge; this is 
easy. Although we do not need the fourth facet, ABG, we include it for the interest of its 
more complicated geometry. 

Of course all the functions Vg and -Rmi depend on whether we are counting cubically or 
affinely (thus, subscripted c or a); the two types will be treated separately. But the general 
conclusions hold that 

R-mi(a;) = 72rs{x) + 36[roAs(a;) + roAcix) + roBcix)] + 12roBix), (32) 

and for the number of reduced symmetry types, rnii(t) with generating function r^\{x). 



rm\{x) = rs(x) + roAsix) + roAcix) + roscix) + rosix), 
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(33) 



where rs(x) is from semimagic and, by Equation @ since |yu(X, y)| = 1 for every lower 
interval in each facet poset (except facet ABC, which we do not use), 

{-l)hoAB{l/x) = Eoab{x) + Eoe{x), 

{-lfroAc{l/x) = EoAcix) + Ead{x) + Ede'{x), ^ ^ 

(34) 

i-lYroBcil/x) = EoBcix) + EoE"ix) + EBoix) + EDE"ix) + Epix), 
{-l)hoB{l/x) = EoBix); 

the sign and reciprocal on the left result from Equation 

There is also the generating function of the number of cubical or affine symmetry classes, 
l{t), whose generating function is l(x). This is obtained from r^\{t) in the same way as L{t) 
is from R^\{t), the exact way depending on whether the count is affine or cubic. 



5.1. Magilatin squares: Cubical count (by upper bound). The weak quasipolynomial 
is exactly as in the semimagic cubical problem. 



5.1.1. Magilatin squares by upper bound. The number of 3 x 3 magilatin squares with strict 
upper bound t is Ldt). We count them via -Rmic(w), the number of reduced magilatin squares 
with largest entry (which we know to be X13) equal to w. The formula is 

t-i 



Lc{t) = J](t - 1 - w)R^Uw). (35) 

w=0 

Equivalently, Rm\c{w) counts ^-integral points in the interior and part of the boundary of 
the inside-out polytope (Qc Jc) of Section WA\ weighted variably by 72/\&s\- 

Now we must calculate the closed Ehrhart generating function for each necessary face. 
This is done by LattE; here are the results. First, OABc'. 

{-ifroAB^^/x) = EoabA^) + EoeA^) 

1 1 

+ 



;i-x)2(l-x2) ^ (l-x)(l-a;)3 (36) 
x + 2 



(1-X)(1-X2)(1-X3) ■ 

Next is OAC^: 

{-lfroAa{l/x) = EoAc^x) + Ead^x) + Ed^x) 

1 1 

+ Tz ^ + 



(1-X)2(1-X2) (1-X2)2 (l_a;2)(l -x3) (37) 

x2 + 2x + 3 

(1-X2)2(1 -x3) ■ 
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The last facet is OBcCc- 

1 1 1 

+ TT^ ^ + 



Finally, the edge 05c: 



2x^ + 1 1 (38) 

(1 — x'^y 1 — 

-2x^ + 4a;^ + 5s + 5 
(1 ■ 

-l)2roB^(l/x) = EoB.(x) = . (39) 



il-xy 

Now Rmic(a;) results from fl32|) . and then from fl35l) we see that 



3v ^ 

Lc(3;) = 7" 77 Rmlc(3;) 

1 — sr 



12x^ 



79x^^ + IQOa;^^ + 260x^^ + 250x^^ + 211x^^ + 179x^° + 181x^ | (40) 
+ 198a;^ + 210a;^ + 181x^ + 125x^ + Qlx^ + 22x^ + 8x^ + Ax + lj 

(l-x4)(l-x5)(l-x3)2(l-x2)2 ' 

The constituents of are extracted as described in Section [?!Tl and here they are: 

{ 3 , 51 4 202 , 3769 , , , . n t ■ 
—t' - —t' + —t' - -—t^ + ci(t)t - Co(t), if t IS even; 
(41) 
3 . 51 4 202 3 7493 , , ^ . ^ t • , , 
— - —t^ + —t' - -^t^ + ci(t)t - co(t), if t IS odd; 

where C\ varies with period 6, given by 

Ci(t) 



'994, 


if t = 


0, 


2 


986, 


if t = 


4 




^ 1909 
2 ' 


if t = 


1 




1925 
2 ' 


if t = 


3, 


5 



(mod 6); 



and Co, given by Table [TJ varies with period 60. Thus the period of turns out to be 60, 
just like that of (not a surprise). However, again as with S'c, the even constant terms have 
half the period of the odd constant terms. That means Lc{2t) has period equal to half the 
denominator of the corresponding inside-out polytope {2P,'K). We have no explanation for 
this. 

The principal constituent, that for t = (mod 60), is 

3 . 51 4 202 3769 n 

—t^ + + 994t - 948. 

10 8 3 10 

The constant term for magilatin squares does not have the simple interpretation as a number 

of linear orderings that it does for magic and semimagic squares, because the entries in the 

square need not all be different. (Still, there is an interpretation as a number of partial 

orderings of the nine cells; see Theorem 4.1 in our general magic and magilatin paper [6], 
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t 


co{t) 


t 


co(t) 


t 


co(t) 


t 


co(t) 


t 


co(t) 





948 


12 


948 


24 


4884 
5 


36 


4452 
5 


48 


4884 
5 


1 


76933 
120 


13 


87301 
120 


25 


16769 
24 


37 


16769 
24 


49 


87301 
120 


2 


2780 
3 


14 


14332 
15 


26 


13036 
15 


38 


14332 
15 


50 


2780 
3 


3 


35607 
40 


15 


6891 
8 


27 


6891 
8 


39 


35607 
40 


51 


32151 
40 


4 


13292 

15 


16 


11996 
15 


28 


13292 
15 


40 


2572 
3 


52 


2572 
3 


5 


18433 
24 


17 


18433 
24 


29 


95621 
120 


41 


85253 
120 


53 


95621 
120 


6 


4452 

5 


18 


4884 
5 


30 


948 


42 


948 


54 


4884 
5 


7 


18497 
24 


19 


95941 
120 


31 


85573 
120 


43 


95941 
120 


55 


18497 
24 


8 


14332 
15 


20 


2780 
3 


32 


2780 
3 


44 


14332 
15 


56 


13036 
15 


9 


32727 
40 


21 


29271 
40 


33 


32727 
40 


45 


6315 
8 


57 


6315 
8 


10 


2572 
3 


22 


2572 
3 


34 


13292 
15 


46 


11996 
15 


58 


13292 
15 


11 


93893 
120 


23 


104261 
120 


35 


20161 
24 


47 


20161 
24 


59 


104261 
120 



Table 7. Constant terms of the constituents of Lc{t), counting all magilatin 
squares by upper bound. 



and recall that an acyclic orientation of a graph can be represented by a partial ordering of 
the vertices.) 

We confirmed the formulas by generating all magilatin squares and comparing the count 
with the coefficients of Lc(x) up to t = 91. 

5.1.2. Symmetry types of magilatin squares, counted by upper bound. The number of sym- 
metry types with strict bound t is /c(^)- We count them via rmiciui), given by Equation 
then 

t-i 

lc{t) = J2{t-l-w)r^,,{w). (42) 

w=0 

Equivalently, r^ic{w) counts the ^-integral points in the interior and part of the boundary 
of the inside-out polytope {Qc, Oc) of Section 14.11 

We get rmic(a;) from ( l33i) : then from ( H2ll we see that 



4 

X 



9x^^ + 20x^^ + 23x^^ + IQx^^ + 10a;" + 13a;^° + 27a;^ + 43a;^l (43) 
+ 54a;'^ + 52x^ + Alx^ + 25x^ + lAx^ + 8x^ + 4a; + 1 
(1 -a;2)2(l -a;3)2(l -x4)(l -x5) 
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The constituents of L are: 



r 1 



Q7 



= { 



240 64 



216 



2029 
720" 



+ ci{t)t — co(t), if t is even; 



^ - — + — - ^l^f + ci(t)t - co(t), if t is odd; 



^.240 64 216 1440 
where Ci varies with period 6, given by 

( 1^ 

2 ' 
151 

144 ' 
131 
V 16 ' 



if t = 


0,2 


if t = 


4 


if t = 


1 


if t = 


3,5 



(mod 6), 



and Co, given by Table [HI varies with period 60. Thus the period of Ic turns out to be 60. As 



t 


co(t) 


t 


co{t) 


t 


co{t) 


t 


co{t) 


t 


co{t) 





9 


12 


9 


24 


47 
5 


36 


41 
5 


48 


47 
5 


1 


49213 
8640 


13 


59581 
8640 


25 


11225 
1728 


37 


11225 
1728 


49 


59581 
8640 


2 


235 
27 


14 


1229 
135 


26 


1067 
135 


38 


1229 
135 


50 


235 
27 


3 


2823 
320 


15 


539 
64 


27 


539 
64 


39 


2823 
320 


51 


2439 
320 


4 


1144 
135 


16 


982 
135 


28 


1144 
135 


40 


218 
27 


52 


218 
27 


5 


12313 
1728 


17 


12313 
1728 


29 


65021 
8640 


41 


54653 
8640 


53 


65021 
8640 


6 


41 
5 


18 


47 
5 


30 


9 


42 


9 


54 


47 
5 


7 


12953 
1728 


19 


68221 
8640 


31 


57853 
8640 


43 


68221 
8640 


55 


12953 
1728 


8 


1229 
135 


20 


235 
27 


32 


235 
27 


44 


1229 
135 


56 


1067 
135 


9 


2503 
320 


21 


2119 
320 


33 


2503 
320 


45 


475 
64 


57 


475 
64 


10 


218 
27 


22 


218 
27 


34 


1144 
135 


46 


982 
135 


58 


1144 
135 


11 


63293 
8640 


23 


73661 
8640 


35 


14041 
1728 


47 


14041 
1728 


59 


73661 
8640 



Table 8. Constant terms of the constituents of lc{t), counting symmetry 
types of magilatin squares by upper bound. 



with Sc and Lc, the period of the even constant terms is half that of the odd constant terms. 
The principal constituent of Ic is 

1 , 3 4 97 o 2029 2 17 

+ + — t - 9. 

240 64 216 720 2 

32 



5.1.3. Some real numbers. For the first several nonzero values of the numbers of magilatin 
squares and of symmetry types, consult this table: 



t 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


Lc{t) 


12 


48 


120 


384 


1068 


2472 


4896 


9072 


15516 


25608 


40296 


61608 


lc{t) 


1 


4 


10 


24 


53 


106 


191 


328 


528 


822 


1230 


1794 


Rm\c{t) 


12 


24 


36 


192 


420 


720 


1020 


1752 


2268 


3648 


4596 


6624 


rmic(i) 


1 


2 


3 


8 


15 


24 


32 


52 


63 


94 


114 


156 



The third line contains the number of symmetry classes of 3 x 3 magilatin squares, counted 
by upper bound. The main numbers, Ldt) and lc(t), are sequences A173548 and A173729 
in the OEIS [T5]. The reduced numbers, -Rmic(^) and r^aicit), are sequences A174018 and 
A174019. In contrast to the semimagic case, the number of squares is not a simple multiple 
of the number of symmetry types. 



5.2. Magilatin squares: AfRne count (by magic sum). The last example is 3 x 3 
magilatin squares, counted affinely. Let -La(t) be the number of 3 x 3 magilatin squares with 
magic sum t > 0. 

The weak quasipolynomial is the same as in affine semimagic. 



5.2.1. Magilatin squares by magic sum. We compute Lg^{t), the number of squares with magic 
sum t, via -Rmia(s), the number of reduced squares with magic sum s. The formula is 

0<s<i-3 
s=t (mods) 

Equivalently, -Rmials) counts ^-integral points in the interior and part of the boundary of the 
inside-out polytope (Qa, ^a) of Section each weighted by 72/\(3s\. 

Now we calculate (by LattE) the closed Ehrhart generating function for each necessary 
face. First, OAB^,: 

{-l)hoABS^/x) = EoabM + Eo£.(x) 

1 1 

~ (l-x)(l-x2)2 + (l-a:)(l-x4) (45) 
2 

~ {1- x){l- x^){l- x^) ■ 

Next is OAC^: 

i-lfroAcA^/x) = EoAcM + EadJo;) + Ed^e'S^) 

_ 1 1 1 

^ (1 -X)(l -X2)2 + (1 -x2)(l-x3) ^ (1 -x3)(l -x4) (46) 

x^ + x"^ + X + 3 

^ (l-x2)(l-x3)(l-x4) ■ 
33 



The last facet is OBgC.^: 



OB^C^ 



x) + EoK"fa;) + E 



E 

+ Eb.s»(x)+E^Jx) 
1 

+ 



1 



(1 - X)(l - X2)2 ^ (1 - x)(l - X4) (1 - x2)(l - X3) 
1 1 

(l-x5)(l-x3)(l-x4)(l + 



1 



X 



Finally, the edge O-Bg 



[-l)2roB,(l/x) = EoBjx) 



1 



(1 -X)(l -X2) 

Now we get Rmia(a;) from (j32|) : then by fj44|) we deduce that 
^3 



X 



1 — x^ 



R-miafa;) 



12x6 < 



' 1 + 3x + Tx^ + ISx^ + 33x^ + 65x^ + 128x6 + 208x^ + 316x^ 
+ 434x^ + 566x^° + 676x" + 784x^^ + 852x^3 + 911x^^ + 936x^^ 
+ 967x^6 ^ 9672;i7 ^ looix^^ + 995x^^ + lOOOx^" + 955x2^ 
^ + 893x^2 + 752x^3 + 624x2^ + 456x^^ + 322x^6 ^ i^^^27 ^ 79^28 

(1 + X)(l + X + X2)(l + x2)(l - x3)(l - x5)(l - x6)(l - x'^){l - X^) 

The constituents of L.^ are the following: 

= -t^ - 3t^ + a2{t)f - ai{t)t + ao{t) - 72Sj{t), 



where 02 varies with period 6, given by 



151 

4 ' 
135 

4 ' 
63 

2 ' 
71 

2 ' 



if t = 
if t = 2, 4 
if t = 1,5 
if t = 3 



(mod 6); 



the linear coefficient varies with period 12, given by 

if t = 



ai(t) 



f 1296 
5 ' 
1347 

10 ' 
831 

5 ' 
2097 

10 ' 
876 

5 ' 
1251 

5 ' 
1257 

10 ' 
2187 
10 ' 



if t = 1,5 
if t = 2, 10 
if t = 3 
if t = 4, 8 
if t = 6 
if t = 7, 11 

if t = 9 

34 



(mod 12); 



the constant term a^, given by Table |9l varies with period 120; and 5*7 is as in the affine 
semimagic count. 



z 


ao[t) 


4- 

Z 


ao[Z) 


Z 


ao[Z) 


Z 


ao[Z) 


4- 

z 


ao[Z) 


Z 


ao[Z) 





876 


20 


340 


40 


400 


60 


840 


80 


376 


100 


364 


1 


4243 
40 


21 


21843 
40 


41 


3283 
40 


61 


5683 
40 


81 


20403 
40 


101 


4723 
40 


2 


1097 
5 


22 


1037 
5 


42 


3597 
5 


62 


917 
~5~ 


82 


1217 
5 


102 


3417 
5 


3 


15219 
40 


23 


—461 
40 


43 


—941 
40 


63 


16659 
40 


83 


— 1901 
40 


103 


499 

lo" 


4 


1604 
5 


24 


4164 
5 


A A 

44 


1484 
5 


r' A 

64 


1784 
5 


O A 

84 


3984 
5 


104 


1664 
5 


5 


1463 
8 


25 


1367 
8 


45 


4887 
8 


65 


1175 
8 


85 


1655 
8 


105 


4599 
8 


6 


3201 
5 


26 


881 
~5~ 


46 


821 
~5~ 


66 


3381 
5 


86 


701 
~5~ 


106 


1001 

5 


7 


3091 
40 


27 


17811 
40 


47 


2131 
40 


67 


1651 
40 


87 


19251 
40 


107 


691 
Iff 


8 


1448 
5 


28 


1388 
5 


48 


3948 
5 


68 


1268 
5 


88 


1568 
5 


108 


3768 
5 


9 


21267 
40 


29 


5587 
40 


49 


5107 
40 


69 


22707 
40 


89 


4147 
40 


109 


6547 
40 


10 


265 


30 


705 


50 


241 


70 


229 


90 


741 


110 


205 


11 


— 1037 
40 


31 


1363 
40 


51 


16083 
40 


71 


403 

lo" 


91 


—77 
'W 


111 


17523 
40 


12 


4092 
5 


32 


1772 
5 


52 


1712 
5 


72 


4272 
5 


92 


1592 
5 


112 


1892 
5 


13 


4819 
40 


33 


19539 
40 


53 


3859 
40 


73 


3379 
40 


93 


20979 
40 


113 


2419 
40 


14 


809 
~5~ 


34 


1109 
5 


54 


3309 
5 


74 


989 
~5~ 


94 


929 
~5~ 


114 


3489 
5 


"1 

15 


4023 
8 


35 


311 
~8~ 


55 


791 
~8~ 


75 


3735 
8 


95 


599 
~8~ 


115 


503 
~8~ 


16 


1676 
5 


36 


3876 
5 


56 


1556 
5 


76 


1496 
5 


96 


4056 
5 


116 


1376 
5 


17 


5011 
40 


37 


7411 
40 


57 


22131 
40 


77 


6451 
40 


97 


5971 
40 


117 


23571 
40 


18 


3273 
5 


38 


593 
5 


58 


893 
5 


78 


3093 
5 


98 


773 
5 


118 


713 
5 


19 


787 
40 


39 


18387 
40 


59 


-173 
40 


79 


2227 
40 


99 


16947 
40 


119 


1267 
40 



Table 9. Constant terms of the truncated constituents of L^{t)^ counting 
magilatin squares by magic sum. 



The period of La turns out to be 840 — the period of the constant terms, due to the 
combination of ao(^) and the constant term of Sj{t). This is equal to the denominator. 
The principal constituent of La, that is, for t = (mod 840), is 

1 4 o 3 151 2 9192 

-t^ - 3t^ + 1 + 948 

8 4 35 

(incorporating the —725*7 term). (As with the cubical magilatin count, there is an interpre- 
tation of the constant term 948 in terms of partial orderings; see Theorem 4.7 in our general 
paper [6].) 

35 



We verified the results by comparing an actual count of magilatin squares with magic sum 
t < 100 to the coefficients in La. 

5.2.2. Symmetry types of magilatin squares, counted by magic sum. We compute la{t), the 
number of symmetry types of squares with magic sum t, via rmia(s), the number of reduced 
symmetry types with magic sum s. The formula is 



Ut) 



E 



if t > 0. 



(51) 



0<s<t-3 
s=t (mods) 



Equivalently, rmia(s) counts ^-integral points in the interior and part of the boundary of the 
inside-out polytope(Qa, Ja) of Section IT2l 

From ( 132|) we get rmia(x) and then from ( 15T1) we see that 

la(x) = -; r^u{x) 



1 

( l + 3x + 7x'^ + I3x^ + 23x^ + 37x^ + 60x^ + 86a;^ + 118a;^ + 149x^'' 
+ 180a;^° + 199x^^ + 212x^^ + 208s^^ + 196a;^^ + 171a;^^ 



^6 + 115x1^ + 96x1® ^ 7g^i9 ^ ^2x^0 ^ q^^2i 



+ QQx^' + 59x^^ + 54x^* + 43x^^ + SSx^*^ + 19x^' + 9x" 

(1 + X)(l + X + X2)(l + x2)(l - x3)(l - x5)(l - x6)(l - X^){1 - X« 

The constituents of /a are: 



(52) 



1^ + a2{t)t^ - ai{t)t + ao{t) - S^{t), 



576 48 

where the quadratic term varies with period 6, given by 

if t = 



02 (t) 



< 



25 
96 ' 
25 

144 ' 
59 

288 ' 



if t = 1,5 
if t = 2, 4 



(mod 6); 



ai{t) 



< 



(mod 12); 



the linear term varies with period 12, given by 

f 31 

15 ' 
103 
120 ' 
133 
120 ' 
47 
30 ' 
37 
30 ' 
233 
120 ' 
11 
15 ' 
203 
120 ' 

and Oo, given by Table [TOl varies with period 120. 

The period of /a is 840. The principal constituent, that for t = (mod 840), is 

1 4 1 . 25 2 74 

1^ 1^ + —t^ 1 + 9. 

576 48 96 35 

36 



if t 


= 




if t 


= 1, 


5 


if t 


= 2, 


10 


if t 


= 3 




if t 


= 4, 


8 


if t 


= 6 




if t 




11 


if t 


= 9 





t 


ao(t) 


t 


do{t) 


t 


doit) 


t 


doit) 


t 


doit) 


t 


doit) 





8 


20 


47 

18 


40 


31 

9 


60 


15 
2 


80 


28 
9 


100 


53 
18 


1 


2027 
2880 


21 


1523 
320 


41 


1067 
2880 


61 


3467 
2880 


81 


1363 
320 


101 


2507 
2880 


2 


553 
360 


22 


493 
360 


42 


257 
40 


62 


373 
360 


82 


673 
360 


102 


237 
40 


3 


979 
320 


23 


-949 
2880 


43 


-1429 
2880 


63 


1139 
320 


83 


-2389 
2880 


103 


11 
2880 


4 


229 
90 


24 


38 
5 


44 


199 
90 


64 


137 
45 


84 


71 
10 


104 


122 
45 


5 


847 
576 


25 


751 
576 


45 


343 
64 


65 


559 
576 


85 


1039 
576 


105 


311 
64 


6 


221 
40 


26 


409 
360 


46 


349 
360 


66 


241 
40 


86 


229 
360 


106 


529 
360 


7 


1739 
2880 


27 


1171 
320 


47 


779 
2880 


67 


299 
2880 


87 


1331 
320 


107 


-661 
2880 


8 


104 
45 


28 


193 
90 


48 


36 
5 


68 


163 
90 


88 


119 
45 


108 


67 
10 


9 


1427 
320 


29 


3083 
2880 


49 


2603 
2880 


69 


1587 
320 


89 


1643 
2880 


109 


4043 
2880 


10 


149 
72 


30 


49 
8 


50 


125 
72 


70 


113 
72 


90 


53 
8 


110 


89 
72 


11 


-1813 
2880 


31 


587 
2880 


51 


1043 
320 


71 


-373 
2880 


91 


-853 
2880 


111 


1203 
320 


12 


73 
10 


32 


131 

45 


52 


247 
90 


72 


39 
5 


92 


217 
90 


112 


146 
45 


13 


2891 
2880 


33 


1299 
320 


53 


1931 
2880 


73 


1451 
2880 


93 


1459 
320 


113 


491 
2880 


14 


301 
360 


34 


601 

360 


54 


229 
40 


74 


481 
360 


94 


421 

360 


114 


249 
40 


15 


279 
64 


35 


-17 
576 


55 


463 
576 


75 


247 
64 


95 


271 
576 


115 


175 

576 


16 


128 
45 


36 


69 
10 


56 


113 
45 


76 


211 
90 


96 


37 
5 


116 


181 
90 


17 


2219 
2880 


37 


4619 
2880 


57 


1491 
320 


77 


3659 
2880 


97 


3179 
2880 


117 


1651 
320 


18 


233 
40 


38 


157 
360 


58 


457 
360 


78 


213 
40 


98 


337 
360 


118 


277 
360 


19 


-277 
2880 


39 


1267 
320 


59 


-1237 
2880 


79 


1163 
2880 


99 


1107 
320 


119 


203 
2880 



Table 10. Constant terms of the truncated constituents of la,{t), the number 
of symmetry types of magilatin squares with given magic sum. 



(This incorporates the —S7 term.) 



5.2.3. Some numbers. The first several nonzero values are given in the table. The third line 
gives the number of symmetry classes of squares. 



t 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


L.it) 


12 


12 


24 


72 


156 


240 


552 


600 


1020 


1548 


2004 


2568 


4008 


4644 


Lit) 


1 


1 


2 


4 


7 


10 


20 


22 


35 


50 


63 


78 


116 


131 


-Rmla(^) 


12 


12 


24 


60 


144 


216 


480 


444 


780 


996 


1404 


1548 


2640 


3696 


Tmlait) 


1 


1 


2 


3 


6 


8 


16 


15 


25 


30 


41 


43 


66 


68 



37 



The sequences Iva(t) and /a(^) are A173549 and A173730 in the OEIS [15j. The numbers 
-Rmia(^) of reduced magilatin squares and rjnia(t) of normahzed, reduced squares with largest 
value t are sequences A174020 and A174021. 

6. Observations and conjectures 

A remarkable fact is that the period of every one of our strong Ehrhart quasipolynomials 
equals the denominator, when it could be much smaller. 

For some small values of t we calculated Ma(t) = Epo j^(t) by hand, which is feasible be- 
cause the problem is 2-dimensional. The process of counting lattice points in a diagram drew 
our attention to some remarkable phenomena that apply to the semimagic and magilatin 
problems as well. Let 6 := dims; let Ck be the coefficient in the quasipolynomial Epo ^(t) 
and let be that in the Ehrhart quasipolynomial of P, and let pk, be their periods. We 
observe that the variation in cs-i is exactly the same as that in c^_i, i.e., 

CS-lit) - - 1) = cj_,{t) - cj_,it - 1), 

but that is not so for most lower coefficients, especially Cq. The reason is that adding each 
new excluded hyperplane results in a constant deduction in degree S — 1 (as we discussed at 
Equation (4.9) in our ffist article |.5j) but a more irregular one in lower terms. We observe 
that Pk increases — that is, there is longer-term variation in Ck — as k decreases in every case. 
Thus we propose some daring conjectures. 

Conjecture 1. In an inside-out counting problem, let S := dim P. 

(a) p"^ \ Pk ioT < k < 6. (We know that ps-i = Pg'^i because cs-i and cj_-^ have the 
same variation.) 

(b) If pj = pj for all j > k, then the variation in Ck is the same as that in c^. 

(c) The period ratios increase by a multiplicative factor as k decreases: 

Pk_ 

Pk 



^ for < < 5. 



Pk-i 

We do not suggest pk \ Pk+i because that is false in general in ordinary Ehrhart theory, 
according to McAllister and Woods |13]. However, it might be true for the kinds of inside- 
out polytopes that arise in cubical and affine counting. 
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